package fem2.material;

import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;

import math2.MutableDouble;
import fem2.Debugger;
import fem2.Element;
import fem2.Variable;
import fem2.element.StructuralElement;
import fem2.enu.EchoLevelType;
import fem2.enu.EquivalentStrainType;
import fem2.enu.State;
import fem2.enu.StressType;
import fem2.enu.VariableName;

/**
 * implicit/explicit version of isotropic damage material model
 * 
 * @author hbui
 * 
 */
public class IsotropicDamageMMimplex extends IsotropicDamageMM {

	protected EquivalentStrainType equivStrainType = EquivalentStrainType.ENERGY_NORM;

	public IsotropicDamageMMimplex(double E, double nu, double rho, double thickness, State ss,
			FluencyCriteria softeningRule, double alpha0) {
		super(E, nu, rho, thickness, ss, softeningRule, alpha0);
	}

	public IsotropicDamageMMimplex(IsotropicDamageMM mm) {
		super(0, 0, mm.getDensity(), mm.getSectionProperty(), mm.getState(), mm.softeningRule,
				mm.alpha0);
		BLAM.copy(mm.De, De);
	}

	/**
	 * set equivalent strain type using in updating the internal variable
	 * 
	 * @param equivStrainType
	 */
	public void setEquivalentStrainType(EquivalentStrainType equivStrainType) {
		this.equivStrainType = equivStrainType;
	}

	@Override
	public IVector getStress(Element e, int k, StressType type) {
		switch (type) {
		case CAUCHY_STRESS:
			return getVariables(e).get(k).getVectorVariable(VariableName.SIGMA);
		case ALGORITHMIC_STRESS:
			return ((StructuralElement) e).getStress(k);
		default:
			throw new Error("invalid stress type");
		}
	}

	@Override
	public void Initialize(Element e) {
		int ng = e.getNumIntegrationPoints();
		int dim = e.getMeshPart().getDimension();
		int ntrs = dim * (dim + 1) / 2;
		ArrayList<Variable> vars = new ArrayList<Variable>();
		for (int i = 0; i < ng; i++) {
			Variable v = new Variable(2);
			v.registerScalarVariable(VariableName.ALPHA, true);
			v.registerScalarVariable(VariableName.STRESS_LIKE_VARIABLE, true);// q
			v.registerScalarVariable(VariableName.DAMAGE, false);
			v.registerVectorVariable(VariableName.SIGMA, ntrs, false);
			v.registerMatrixVariable(VariableName.TANGENT_OPERATOR, ntrs, false);
			vars.add(v);
		}
		setVariables(e, vars);
		/*
		 * Initialize constitutive tangent for each element at each Gauss point
		 */
		for (Variable v : vars) {
			/*
			 * initialize tangent operator
			 */
			IMatrix C = v.getMatrixVariable(VariableName.TANGENT_OPERATOR);
			BLAM.copy(this.De, C);
			/*
			 * initialize internal variable and stress-like variable
			 */
			MutableDouble alpha = v.getScalarVariable(VariableName.ALPHA);
			MutableDouble q = v.getScalarVariable(VariableName.STRESS_LIKE_VARIABLE);
			alpha.setValue(alpha0);
			q.setValue(alpha0);
			/*
			 * save the variable to create history
			 */
			v.save();
			v.save();
		}
	}

	@Override
	public void InitializeSolutionStep(Element e) {
		StructuralElement se = (StructuralElement) e;
		for (int k = 0; k < se.getNumIntegrationPoints(); k++) {
			double[] xi = se.getIntegrationPoint(k);
			IVector sigma = se.getStress(k);
			IMatrix C = getConstitutiveMatrix(e, k);
			IVector epsilon = se.computeStrain(xi);

			// if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
			// Debugger.watch("initialize step:");
			// Debugger.watch("epsilon =" + epsilon);
			// Debugger.watch("sigma =" + sigma);
			// }

			/*
			 * perform first explicit update by computing extrapolated internal
			 * variable
			 */
			double alpha_old_old = getVariables(e).get(k).getScalarVariable(VariableName.ALPHA, -2);
			double alpha_old = getVariables(e).get(k).getScalarVariable(VariableName.ALPHA, -1);
			double dalpha = alpha_old - alpha_old_old;
			double alpha_alg = alpha_old + dalpha;

			double q_old = getVariables(e).get(k).getScalarVariable(
					VariableName.STRESS_LIKE_VARIABLE, -1);
			double H = softeningRule.computeFluencyDerivative(alpha_old);
			double q_alg = q_old + H * dalpha;
			/*
			 * compute damage based on extrapolated value
			 */
			double damage_alg = 1 - q_alg / alpha_alg;
			/*
			 * compute extrapolated stress
			 */
			BLAM.multiply(1 - damage_alg, BLAM.NO_TRANSPOSE, De, epsilon, 0.0, sigma);
			/*
			 * compute extrapolated algorithmic tangent
			 */
			BLAM.copy(De, C);
			BLAM.scale(1 - damage_alg, C);
		}
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		StructuralElement se = (StructuralElement) e;
		int ng = e.getNumIntegrationPoints();

		for (int k = 0; k < ng; k++) {
			double[] xi = e.getIntegrationPoint(k);
			IVector sigma = se.getStress(k);
			IVector corrected_sigma = getVariables(e).get(k).getVectorVariable(VariableName.SIGMA);
			IMatrix C = getConstitutiveMatrix(e, k);
			MutableDouble alpha = getVariables(e).get(k).getScalarVariable(VariableName.ALPHA);
			double alpha_old = getVariables(e).get(k).getScalarVariable(VariableName.ALPHA, -1);
			double alpha_old_old = getVariables(e).get(k).getScalarVariable(VariableName.ALPHA, -2);
			MutableDouble q = getVariables(e).get(k).getScalarVariable(
					VariableName.STRESS_LIKE_VARIABLE);
			double q_old = getVariables(e).get(k).getScalarVariable(
					VariableName.STRESS_LIKE_VARIABLE, -1);
			MutableDouble damage = getVariables(e).get(k).getScalarVariable(VariableName.DAMAGE);

			/**********************************************
			 * perform implicit update to compute corrected stress & damage
			 **********************************************/
			/*
			 * compute strain
			 */
			IVector epsilon = se.computeStrain(xi);

			/*
			 * compute effective stress
			 */
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, De, epsilon, 0.0, sigma);

			/*
			 * compute equivalent strain
			 */
			double etilde;
			switch (equivStrainType) {
			case ENERGY_NORM:
				etilde = Math.sqrt(BLAM.dot(epsilon, sigma));
				break;
			default:
				throw new Error("invalid equivalent strain type");
			}

			if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
				Debugger.watch("etilde = ", etilde);
			}

			/*
			 * compute trial state
			 */
			double alpha_trial = alpha_old;

			/*
			 * check loading/unloading criteria (use g to be consistent with the
			 * paper)
			 */
			double g = etilde - alpha_trial;
			double dlambda;
			if (g <= 0.0) {
				/*
				 * unloading/elastic state
				 */
				dlambda = 0.0;
			} else {
				/*
				 * loading/damage state
				 */
				dlambda = g;
			}

			/*
			 * update internal variable
			 */
			alpha.setValue(alpha_old + dlambda);

			/*
			 * update stress-like variable
			 */
			double H = softeningRule.computeFluencyDerivative(alpha_old);
			// double H =
			// softeningRule.computeFluencyDerivative(alpha.getValue());
			q.setValue(q_old + H * dlambda);

			/*
			 * update damage variable
			 */
			damage.setValue(1 - q.getValue() / alpha.getValue());
			// if (damage.getValue() > 0.9999) {
			// damage.setValue(0.9999);
			// }
			if (damage.getValue() < 0.0) {
				/*
				 * just for debugging, to make sure the input & material model
				 * work correcly
				 */
				throw new Error("negative damage. It's unphysical");
			}

			/*
			 * update stress
			 */
			BLAM.copy(sigma, corrected_sigma);
			BLAM.scale(1 - damage.getValue(), corrected_sigma);

			/**********************************************
			 * perform explicit update to compute algorithmic stress and tangent
			 **********************************************/

			double dalpha = alpha_old - alpha_old_old;
			double alpha_alg = alpha_old + dalpha;

			// H = softeningRule.computeFluencyDerivative(alpha_old);
			double q_alg = q_old + H * dalpha;
			/*
			 * compute damage based on extrapolated value
			 */
			double damage_alg = 1 - q_alg / alpha_alg;
			/*
			 * compute extrapolated stress
			 */
			BLAM.scale(1 - damage_alg, sigma);
			/*
			 * compute extrapolated algorithmic tangent
			 */
			BLAM.copy(De, C);
			BLAM.scale(1 - damage_alg, C);

			/**********************************************
			 * debugging part
			 **********************************************/
			if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
				Debugger.watch("finalize step:");
				Debugger.watch("epsilon =" + epsilon);
				Debugger.watch("sigma =" + sigma);
				Debugger.watch("alpha = ",
						getVariables(e).get(k).getScalarVariable(VariableName.ALPHA));
				Debugger.watch("damage = ",
						getVariables(e).get(k).getScalarVariable(VariableName.DAMAGE));
			}
		}
	}

}
